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Oh. We develop a formalism for the calculation of the frequency band structure of a 

■ phononic crystal consisting of non-overlapping elastic spheres, characterized by Lame 

coefficients which may be complex and frequency dependent, arranged periodically 
in a host medium with different mass density and Lame coefficients. We view the 
crystal as a sequence of planes of spheres, parallel to and having the two dimensional 
periodicity of a given crystallographic plane, and obtain the complex band structure of 
the infinite crystal associated with this plane. The method allows one to calculate, also, 
C/3 ' the transmission, reflection, and absorption coefficients for an elastic wave (longitudinal 

or transverse) incident, at any angle, on a slab of the crystal of finite thickness. We 
demonstrate the efficiency of the method by applying it to a specific example. 



43.20.+g, 43.40.+S, 46.40.Cd, 63.30.+d 



'. ■ I. INTRODUCTION 

The elastic properties of a locally homogeneous and isotropic composite material are characterized by a mass density 
O | p and Lame coefficients A and fj, which vary in space.u The composite materials, we shall be concerned with in this 
paper consist of homogeneous particles (solid or fluid inclusions the dimensions of which must be large enough in 
| order for a macroscopic description of their elastic properties to be valid) distributed periodically in a host medium 
. characterized by different mass density and Lame coefficients. We assume, throughout this paper, that the particles 
do not overlap with each other (cermet topologyo). The alternative case, when the particles connect with each other 
to form a continuous network is also interesting but will not concern us here. When identical particles are distributed 
£0 ' periodically in a host medium, the composite material may be referred to as a phononic crystal. In this case the mass 
^■j^ , density and the Lame coefficients vary periodically in space: 

p(r + R„) =p(r), n(T + K n )=n(T), A(r + R„) = A(r), (1.1) 

■ where {R„} denotes a Bravais lattice. 

In recent years there has been a growing interest, in the study of phononic crystals which is inspired to a large 
degree by corresponding work in photonic crystals Bu These are composite materials with a dielectric function which 
varies periodically in space. A typical example: identical particles, large enough to be describable by a macroscopic 
dielectric function, are arranged periodically in a host material with a different dielectric function. Photonic crystals 
have many interesting properties both in relation to basic physics and technological applications. In particular, the 
existence of absolute frequency gaps (photonic gaps) in certain such crystals, i.e., regions of frequency over which 
electromagnetic (EM) waves can not exist within the crystal, has attracted a lot_pf attention, mainly because of 
promising applications in optoelectronics, as pointed out initially by Yablonovich.B In principle, one can design a 
perfect mirror, non-absorbing over a selected region of frequency (corresponding to a photonic gap), a non-absorbing 
resonance cavity, etc.Jj A number of theoretical calculations, predict the existence of such gaps for appropriately 
designed photonic crystals, but so far only crystals which exhibit gaps up to the infrared region have been constructed.Q 
However, progress to higher frequencies is expected in the near future. In relation to basic physics, photonic crystals 
are interesting in a number of ways.H For example, they can be the starting point hipa process of gradual introduction 
of disorder and a study of consequent phenomena, including Anderson localization.!! 

—New, phononic crystals have properties which mirror those of photonic crystals and corresponding applications too 
.t2rta With an appropriate choice of the parameters involved one may obtain phononic crystals with absolute frequency 
gaps (phononic gaps) in selected regions of frequency. An elastic wave, whose frequency lies within an absolute gap of 
a phononic crystal, will be completely reflected by it; from which follows the possibility of constructing non-absorbing 
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mirrors of elastic waves and vibration-free cavities which might be very useful in high-precision mechanical systems 
operating in a given frequency range. A*id in relation to basic physics, one can use elastic waves to study phenomena 
such as those associated with disorder,E3 in more or less the same manner as with EM waves. 

There are, however, some essential differences between EM and elastic waves and this means that the normal modes 
of the elastic field in a phononic crystal are in some ways quite different from those of the EM field in a photonic 
crystal. In a homogeneous isotropic medium the elastic waves can, in general, be purely longitudinal (in which case 
the displacement vector u(r) satisfies the condition V x u = 0) or purely transverse (in which case V • u = 0). In a 
phononic crystal this is no longer the case and a normal mode usually has a longitudinal and a transverse component. 
One expects that because of this coupling between longitudinal and transverse waves, it will be more difficult to obtain 
absolute frequency gaps in a phononic crystal. We recall that the normal modes of the EM field in a photonic crystal 
are exclusively transverse. On the other hand there are, in general, more parameters relevant to the determination 
of phononic gaps than there are in the determination of photonic gaps. In the case of a binary system (consisting 
of material 2 distributed in material 1) we have for photonic crystals two independent parameters: the ratio of the 
dielectric functions of the two materials £2/61; and the fractional volume occupied by material 2, to be denoted by 
/. For phononic crystals there are five independent parameters: 112/1*1, A2/A1, pi/ pi, and /; where pj, Xj, jij 

denote the mass density, and the Lame coefficients of material j — 1,2. 

We have, so far, implicitly assumed that the Lame coefficients describing the constituent materials of the phononic 
crystals are all different from zero, real quantities, and constant (independent of the frequency). But this is not always 
the case. The phononic crystal may consist, for example, of solid particles (material 2) arranged periodically (at least 
approximately) in a liquid (material 1). If the liquid is a normal fluid, like water, fix — and the transverse sound 
in the liquid is suppressed. This, however, is not the case, for a viscous fluid. The role of shear viscosity in phononic 
crystals has been pointed out by Sprik and Wegdam. Shear viscosity is equally important in phononic crystals 
consisting of liquid particles in a solid host background (liquid-containing porous solidsEira). Colloidal suspensions 
of solid spheres in a liquid, also, have interesting acoustic properties. cB Finally, it may be of some interest to consider 
composite materials consisting of two liquids (e.g. drops of oil in water) although in this case a periodic arrangement 
of the drops can only be a rough approximation to the real system. Ji appears that acoustic gaps are easily obtained 
in three-dimensional (3D) fluid-fluid composites, when [i\ = p,2 = 0.cil 

The few calculations published so far relating to 3D phononic crystals deal, almost exclusively, with the fp^qucncy 
band structure of these crystals, which is obtained via a plane-wave expansion of the displacement field.E3lij On 
the other hand, a lot of theoretical and experimental work has been done on systems with two-dimensional (2D) 
periodicity, with translational invariance along the third dimension. A typical example of such systems consists of 
a set of long identical cylinders parallel to the z direction, crossing the xy plane at the sites of a-2D-Jattice. By 
considering waves propagating normal to the cylinders, the problem is reduced to two dimensions. L3t2l The above 
investigations have shown that phononic gaps are possible in both 2D and 3D systems. 

Although knowing the frequency band structure of a phononic crystal is very useful, more is required for a full 
interpretation and analysis of the experimental data. In an experiment one usually measures the reflection and/or 
transmission coefficients of an acoustic / elastic wave incident on a slab of the phononic crystal, and consequently theory 
should be able to provide reliable estimates of these, the experimentally measured quantities, as well. The so-called 
on-shell methods developed injelation to photonic crystals can do exactly that, besides an accurate evaluation of 
the frequency band structure. E3~E3 In these methods one determines for a given frequency u> and a given reduced 
wavevector, k|| , parallel to a given crystallographic plane of the crystal, the Bloch-wave solutions of the elastic field of 
the infinite crystal; these consist of propagating and evanescent waves. The propagating waves constitute the normal 
modes of the infinite phononic crystal. The evanescent waves do not represent real waves, they are mathematical 
entities which enter directly or indirectly (depending on the method of calculation) into the evaluation of the reflection 
and transmission coefficients of a wave, with given lo and k|| , incident on a slab of the crystal parallel to the given 
crystallographic plane. On-shell methods have certain advantages over the plane-wave method, even if one is only 
interested in the frequency band structure and the corresponding normal modes of vibration of the infinite phononic 
crystal. In an on-shell method one can easily allow the Lame coefficients of any of the constituent materials of the 
crystal to depend on the frequency, as is necessary in some cases, without any difficulty, which is not the case with 
the plane-wave method. And, as a rule, on-shell methods are computationally more efficient .c3'E3 

The on-shell method we describe in the present paper is analogous to that which some of us have developed for 
photonic crystals. LJ It applies to systems which consist of non-overlapping spherical particles arranged periodically 
in a host medium characterized-.by different ma ss density and Lame coefficients. Sections [n] to VI are devoted to the 
development of the formalism.!^] In section VII we demonstrate the applicability of the method on a specific system: 
an fee crystal of silica spheres in ice. Finally the last section concludes this article. 
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II. MULTIPOLE EXPANSION OF THE ELASTIC FIELD 



The displacement wector U(r, t), in a homogeneous elastic medium of mass density p and Lame coefficients A, /i, 
satisfies the equation^ 



(A + 2/i) V(V-U)-/iVx(VxU)- pd 2 t V = 0. 
In the case of a harmonic elastic wave of angular frequency lo, we have 

U(r, *) = Re [u(r)exp(-iurf)] , 



and Eq. (2.1) reduces to the following time- independent form 

(A + 2/i) V (V • u) - (iV x (V x u) + puj 2 u = 0. 



(2.1) 
(2.2) 
(2.3) 



We note that for ordinary elastic media the Lame coefficients are real numbers. Media,where loss is possible, assuming 
the time dependence given in Eq. (|2.2[), are described by complex Lame coefficientsLil: 



A = A e — iujXy, n = fi e — \wfi v . 



(2.4) 



The most general solution of Eq. (2.3) consists of two elastic waves which propagate independently. These are: a 
longitudinal (irrotational) wave, which satisfies the equations 



V 2 u + qf u = 0, V x u = 0, 



(2.5) 



where qi = uj/ci, q = + 2/i)/p being the speed of propagation of this wave; and a transverse (divergenceless) 
wave, which satisfies the equations 



V 2 u + q 2 t u = 0, V • u = 0, 
where qt = w/ct, Ct = \J p/p being the speed of propagation of this wave. 



(2.6) 



In the present paper we shall use, besides t he m ore familiar solutions of Eqs. (2.5) and (p.6|) representing longitudinal 



and transverse plane elastic waves (see Eq. ( J3. If) be low), the so-called spherical- wave solutions of these equations. A 
complete set of spherical- wave solutions of Eq. (2.5), known as irrotational vector wave functions, is given byE3 



u L( r ) 



-V[//(«nr)17*(f)] 



(2.7a) 



where fi may be any linear combination of the spherical Bessel function, jg, and the spherical Hankcl function, lit. 
Y™(r) are the usual spherical harmonics, with f denoting the angular variables (0, <f>) of r in a system of spherical 
coordinates. I— . 
A complete set of spherical- wave solutions of Eq. (|2.6|) is given byL 2 ! 



nf m (r) = Mq t r)X em (r) 



and 



O) = -V x h{q t r)yUm{i) , 



(2.7b) 



(2.7c) 



which are also known as solenoidal vector wave functions. The vector spherical harmonics, denoted by X£ m (f), are 
defined by 



y/e(£+l)Xt m (r) = L!7 l (f) = -ir x VY e m (r) 
By definition Xoo(f) = 0; for £ > 1 we have 



y/e(e+l)Xe m (r) = [a- m cos0 e 1 * Y™- 1 ^) -msind Y e m (r) + af cos6» e _1 * Y e m+1 (r)] e e 
+ i [aj m Y e m -\r) - af Y £ m+1 (r)] e , 



(2.8a) 



(2.8b) 



where 
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aT = ^[{e~m)(£ + m+l)} 1/2 , (2.8c) 

and eg, e^, are the usual polar and azimuthal unit vectors, respectively, in the chosen system of spherical coor di nates . 

The most general displacement field can be written as a linear sum of the spherical waves given byEqs. (^7a|)-(|2~7c1), 
as follows 

u(r) = V (af„/,(g t r)X £m (f) + a^-V x Mq t r)X im (r) + a^-V [Mqir)Y e m (r )}} , (2.9) 
where a»L, P = M, iV, L, are coefficients to be determined. 

III. SCATTERING OF A PLANE WAVE BY A SPHERE 

A plane elastic wave, of wavevector q, propagating in a homogeneous elastic medium has the form 

Ui n (r) = uo(q) exp(iq • r), (3.1) 

with u (q) = Mo(q)e, where u Q denotes the magnitude and e, a unit vector, the polarization of the displacement field. 
In the case of a longitudinal plane wave we can write: q = qie q and e = e q . Since the plane wave is finite everywhere, 



its multipole expansion into spherical waves, according to Eq. (2.7a), involves only the radial functions je(qir); we 
have 



u in (r) = ]T a^-V \jt(qir)Yt(r)] . (3.2) 

One can easily show that the coefficients a^ are given by 

ag£ = A$£(q).uo(q), (3.3) 

where 

A^(q) = 4 7 ri^ 1 (-l) m+1 y £ -™(q) e g . (3.4) 
In the case of a transverse plane wave w e ha ve: q = qte q and e_Le 9 . Such a wave can be written as a linear sum of 



the spherical waves given by Eqs. fl2.7b| ) and ( 2.7c ), and again involves only the radial functions jt(qtr); we have 



M r ) = J2{^Mqtr)X em (r) + a^V x u{q t r)X tm {v)\ . (3.5) 

The coefficients ai?, with P = M, N, can be written as 

a^=A^(q)-Uo(q), (3.6) 



where 



and 



x | [of cos 6» e 1 * ^- m - 1 (q) + msin(9 Y e " m (q) + o^™ cos0 e _i * *7 m+1 (q)] e e 

+i [a? Y- m -\q) - a- m e"* Yf m+1 (q)] e^j , (3.7) 

A \m+l ( 

- [aj» cos<9 y / _m-1 (q) + m sin F^ m (q) + V" cos<9 e'^ Y^ m+1 (q)] e | , (3.8) 
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where 9 and <j) denote the angular variables of q in the chosen system of spherical coordinates. 

We now consider a sphere of radius S, centered at the origin of coordinates. We assume that the sphere, which has 
a uniform mass density p s , is embedded in a homogeneous medium of mass density p; the wavenumbers of the elastic 
waves in the sphere (q sv ) and in the host medium (q u ), where v = I or t, are also different. When a plane wave is 
incident on the sphere, it is scattered by it, so that the wavefield outside the sphere consists of the incident wave a nd a 
scattered wave. Since the scattered wave is outgoing at infinity, its expansion in spherical waves is given by Eq. (|2.9| ) 
with fi = ht, which has the asymptotic form appropriate to an outgoing spherical wave: h~f(qr) » (— i) e exp(iqr)/iqr 
as r — ► oo. We have 



u sc (r) = V(a+f h+(q t r)X lm (r) + a+^-V x hj{q t r)X lm {v ) + a+^-V [hj ( qi r)Y e m (r ) 



(3.9) 



The wavefield inside the sphere is given by Eq. (2.E) with fp = ji, since it must be finite at the origin; we have 

^(a™j £ (fer)X fa (r) + a^— V x n {q st r)X (m {v) + a&— V \jt(q s ir)Y e m (r)]l . (3.10) 



The coefficients a^, a^, P = M, N, L, in Eqs. ( |3.9| ) and ( [3.10 ) are determined by the requirement of continuity of 
the displacement vector, u(r), and of the surface traction, r(r) = <r(r) • f , at the surface of the sphere; <r(r) denotes 
the stress tensor. The components of the surface traction are given by (see, e.g., Ref. fy) 



r r = AV • u + 2pd r u r 
re = p 

T 

1 



-00U r + d r v,0 

r r 



t~4> = M 



r sin 6 



OrU^ ^ 



(3.11a) 
(3.11b) 

(3.11c) 



The continuity of u r , ug, w^, r r , Tg, at the surface of the sphere allows us to determine uniquely the coefficients 
& tm C = M, N, L) of the scattered wave, given by Eq. (3.9), and the coefficients a|^ of the wa ve in side the sphere, 
given by Eq. ( [3.10 ), in terms of the known coefficients eq^ of the incident wave, given by Eqs. (3^) or (|3~5|). After 
some lengthy but straightforward algebra one obtains (see, e.g., Ref. [H]) 



PP' OP 
'm:£'m' ^I'm 



(3.12) 



Explicit expressions for the nonzero elements of the T matrix in the case of a solid scatterer in a solid host are given 
in Appendix Similar expressions for the cases involving a liquid scatterer or host can be found in Ref. [5l]. 



IV. SCATTERING BY A PLANE OF SPHERES 



We consider a plane of spheres at z = 0: an array of spheres centered on the sites of a 2D lattice specified by 

R„ = rtiai + n 2 a 2 , (4.1) 



where ai and a 2 are primitive vectors in the xy plane and rii, ti 2 = 0, ±1, +3, . . .. 
The corresponding 2D reciprocal lattice is obtained in the usual manneioO as follows 

g = ?nibi + m 2 b 2 , 

where mi, m 2 = 0, ±1, ±2, ±3, . . . and bi, b 2 are defined by 

b l • slj = 2nSij . 



(4.2) 



(4.3) 



We now assume that a plane wave (it can be longitudinal or transverse) is incident on the plane of spheres. We 
write the displacement vector u; n (r) corresponding to it as follows 



u ta( r ) = E t Ui »]gV expCiKg'*,' • r)e^ , 



(4.4) 
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where s' = +(— ) corresponds to a wave incident on the plane of spheres from the left (right); v 1 specifies the 
polarization of the incident wave: q v i = qi = u/ci for a longitudinal wave and q v > = q t = Lu/c t for a transverse wave; 



K 



^k i|+ g'±[e-( k||+ g') 2 ] 1/2 e 



g v 



(4.5) 



where e z is the unit vector along the z axis, and we have written the component of the incident wavevector parallel 
to the plane of spheres as the sum of a reduced wavevector kn, which lies in the surface Brillouin zone (SBZ) of the 
given lattice, and an appropriate reciprocal-lattice vector g'. This is always possible and it facilitates the subsequent 
calculation. For v' = I, i' = 1 denotes the only non-zero component of the displacement vector, ei being the radial 
unit vector along the direction of K g/; . For v' = t, i' = 2, 3 denote the only non-zero components of the displacement 

vector; e2,e3 being the polar and azimuthal unit vectors, respectively, which are perpendicular to K g , t . In the same 



manner (as in Eq. (|4.5|)) we define, for given ky, uj, a wavevector K g „ and the corresponding for any g and any v 



We remember that the i = 1 component of the displacement vector is always associated with a longitudinal plane wave 
(y = I), and that the i = 2,3 components of the displacement vector are always associated with a transverse plane 
wave (y — t), so that the character (longitudinal or transverse) of a given plane wave is automatically determined 
by the non- vanishing components of the displacement vector associated with it, and need not be stated explicitly in 
every case. When (ky + g) 2 > q 2 the corresponding wave decays to the right for s = + , and to the left for s = — ; and 
the corresponding unit vectors become complex. Indeed, the unit vectors arc defined in a Cartesian system of 
coordinates as follows 



ei 



e T sm u cos ( 



e„ sm f sm ( 



where and cj> denote the angular variables of K g/ , and 



^2 = e-r COS t> COS (f> - 

£3 = — e x sin <p + e v 



1 sm < 



e 2 sm ( 



(4.6a) 



(4.6b) 
(4.6c) 



where 8 and </> here denote the angular variables of K gt . We note that the z component of K gt/ (denoted by K^ vz ) is 



real if (k|| -l-g) < qf, and imaginary if (kii 



-g) 2 > <&■ 



In the latter case, cos#k s in Eqs. (4.6a) and (4.6b) is replaced 



by K s %vz jq ll and sin^Kj,, by |k|| + g\/q v , so that ei and 62 become complex. 

Beca use o f the 2D periodicity of the structure under consideration, the wave scattered from it, when the wave given 
by Eq. (^4) is incident upon it, has the form 



u sc (r) = b^ J ^exp(ik|| • R n )hj (q t r n )X im (in) 

ira R„ 

+ h em— V x y^ ex P( ik || ' R n )hj (q t r n )X em (r n ) 
+ b £ + ™- V 5> x P( ik H • Rn)hj( qi r n )Y e m (r, 



R„ 



(4.7) 



where r„ = r — R„. We note that exp[i(k|| + g) • R„] = exp(ikii • R„) because of Eq. (O). Eq. (4/7) tells us that the 
scattered wave is a sum of outgoing spherical waves centered on the spheres of the plane, and that the wave scattered 
from the sphere at R„ differs from that scattered from the sphere at the origin (R„ = 0) only by the phase factor 
exp(ik|| • R n ). We note the presence in the scattered wavefield of both longitudinal and transverse waves even when 
the incident wave is purely longitudinal or purely transverse. 

The coefficients which determine the scattered wave from the sphere at the origin are determined from the 
total incident wave on that sphere, which consists of the incident plane wave and the sum of the waves scattered from 
all the other spheres in the plan e. T he latter, denoted by Ug C (r), is obtained from u sc (r) by the removal of the term 
corresponding to R„ = in Eq.(4/7). Ug C (r) can be expanded into spherical waves about the origin as follows 



,(') = E 



(in 



(f)+b^-Vxj / ( 9t r)X <m (f) 
Qt 



qi 



It can be shown (see Appendix ph that 



= E « 

P'l'm 1 



ppi 



(4.8) 



(4.9) 
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It is worth noting that the matrix elements of ft depend on the geometry (4.1) of the plane and, through q v , on 
the frequency, the mass density and the Lame coefficients of the medium surrounding the spheres of the plane; they 
depend also on the reduced wavevector k|| of the incident wave; but they do not depend on the scattering properties 
of the individual sphere. 

The coefficients b^, which describe the scattered wave from the sphere at the origin of the coordinates, are given 

by 



b +p - 

u im — 



P't'm 



T 

J- 1 



PP 

tm;t 



I'm' 1 A t'm' T - Ugr m i 



(4.10) 



The coefficients on the right-hand side of Eq. (4.10) describe the total wave incident on the sphere at the origin of 



coordinates; a^ derive from the incident plane wave given by Eq. Q via Eqs. Q and Q, and b p n from the 
field defined by Eq. (Q. Combining Eqs. (fl9|) and Q4.10| ), we obtain 



E 

P'l'm' L 



5ppi 



_ \p t pp" Q p"p' 

/ j J -im;l"m" ^ L I" m" ;l' m' 



P"£"m" 



b& = E 

P'l'm 1 



T 



PP' 



OP' 



(4.11) 



Eq. (4.11) determines the coefficients bt of the wave scattered fr om the pl ane o f spheres, given by Eq. (4.7), in terms 



of the coefficients a^ of the inciden t w ave. According to Eqs. (^3) and (3J3), we write the coefficients a"^ of the 
incident plane wave, defined by Eq. (4.4), in the form 



Jg'i' ; 



(4.12) 



where A° p are given by Eqs. (3.4), (3.7), and fl3.8|). Due to the linearity of Eqs. (4.11), the coefficients b^ p can be 
written as follows 



— E ^<m;i'(^gv) [ u ™]gV > 



(4.13) 



so that the system of Eqs. ( 4.11 ) reduces to 



E 



P'l'm' . 



SpP'Su'S-n 



V T PP " Q P " P ' 

/ j ± lm;l"rn" iL l"m"\t'm' 



o+P' /fra' \ _ rpPP' ,0P' iffs' 



(4.14) 



We remember that i' , s' , and g' are parameter values characterizing the incident wave (we remember also that v 1 is 
determined by i': v' = / for %' = 1 and v' = t for i' = 2,3). Eqs. (4.14) (or, equivalently, Eqs. (4.11)) constitute a 
system of infinitely many linear equations. It is solved by introducing an angular momentum cut-off, ^ max , truncating 
all angular momentum expansions to £ max , thus reducing the dimension of the system to 3^ ax + 6£ max -I- 1. Moreover, 
by using the properties of the matrix elements Qfm-t'm' gi ven by Eqs. ( Bll ) of Appendix this system can be 
reduced to two independent systems of (3ljt av + 5£ max )/2 and (3£^ ax + 7£ max + 2)/2 linear equations, respectively. 
Finally, the scattered wave given by Eq. (4.7) can be expressed as a sum of plane waves using the following identity 



exp(ikii • R n )hj(q^r n )Y e m (i ri 



2 ' T( -y/"(K±)exp(iK± -r) 



(4.15) 



where Aq denotes the area of the unit cell of the lattice given by Eq. ([O]). The plus (minus) sign on K gl/ must be 
used for z > (z < 0). We note that K^ vz can be real or imaginary. In the latter case cos# K ± in the standard 



formulae for Y™(K .£„) is replaced by K^ vz jq v (see text following Eq. ( 4.6c[ )). 

Using Eq. (4.15) we can expand the scattered wave into a series of longitudinal and transverse plane waves, as 
follows 



<( r ) = E E b ^T A L(K g J exp(iK g „ ■ r) , 

g Plm 



(4.16) 



where 
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(4.17a) 



27r(-i)< 



c| [«7 m costf e 1 * F™- 1 ^) - msintf F/"(K gt ) + a? cosd e" 1 * Y t m+1 (K s st )] e 2 
-i [a t e v Y t (K gt ) - o £ e ^ (K gt )J e 3 J. , 



(4.17b) 



27f(-i)' 



= h[or t m j*Y™-\Klt) 



aT e -i* y-+i(K|J] e 2 



[«7 m cos^ F/"- 1 ^) - msinfl F/^K^) + "7* cosfl e" 1 * y™ +1 (K gt )] e 3 J 



(4.17c) 



with 9 and denoting the angular variables of K gt . Substituting h~J m from Eq. (4.13) into Eq. (4.16) we obtain 



u fc( r ) = E t Wsc ]g* ex P( iK g^ • r ) e» , 



g< 



with 



j g , 



= E E ^(KJJB^K^,) [u m ] g ' v , 



(4.18) 



(4.19) 



i' P£m 



where the superscript s = +(— ) holds for z > (z < 0). We note that the K gJ/ in Eq. (|4.18| ) have the same frequency 



u> (the same wavenumber q u ) and the same reduced wavevector k|| as the incident wave. We remember that for i = 1 
Z and, for a given g, ei is the radial unit vector along the direction of K g( . Similarly, for i = 2, 3, v — t and, for 



given g, e2,e3 are the polar and azimuthal unit vectors, respectively, which are orthogonal to K gt . Eq. (4.18) tells 
us that the scattered wave consists, in general, of a number of diffracted beams, of the same w and k|| , corresponding 
to different g vectors and polarization modes (longitudinal or transverse) . We note, however, that only beams for 
which Kg VZ is real constitute propagating waves. The coefficients in Eq. ( 4.18| ), given by Eq. ( 4.19 ), are functions 



of the Bg^.j, (K gV ,) coefficients and through them depend on the incident plane wave. These coefficients are to be 
evaluated for an incident longitudinal (i' = 1) or transverse (%' = 2, 3) plane wave, with a wavevector K g/J/ , given 



by Eq. (4.5), incident from the left (s 1 = +) or from the right (s' = — ), with a displacement vector along the i'th 



direction of magnitude equal to unity. In other words, B^ i( (K^) are calculated from Eq. (4.14), substituting 



4w(«w), given by Eqs. @, ([[]), and (g 



on the right-hand side of this equation. Obviously, when i 1 = 1, 

only the coefficients ^.^(K*^) are nonzero, and when i' = 2, 3 only -^m-i'^s't) an( l ^Sn-,v C^g't) are nonzero. 

Let us for the sake of clarity assume that a plane wave given by Eq. ( |4.4| ) is incident on the plane of spheres from 
the left as in Fig. [pa. Then the transmitted wave (incident+scattered) on the right of the plane of spheres can be 
written as 



u t + r( r ) = E ex P( iK ^ ' r )^ ' z > °> 

gi 



with 



Kr] gi = [ttfajji <W + Kc] gi = E M gii'i' 



U in] g 'i' 



(4.20) 



(4.21) 



and the reflected wave as 



u rf ( r ) = E [ Urf ]gi ex P( iK g» ' > z < °' 



(4.22) 



g< 



with 
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M gi = [u sc ] sl = M S ifg'i> faiJg'i' 



(4.23) 



Eqs. ([4.19[) , ( J4.2l|) and (|4.23|) define the transmission (M++) and reflection (M _+ ) matrix elements for a plane wave 
incident on the plane of spheres from the left. Similarly, we can define the transmission matrix elements M^^,^ and 

the reflection matrix elements M~^7 &li , for a plane wave incident on the plane of spheres from the right (see Fig. ^b). 
We obtain 



(4.24) 



Ptm 



One can show that the above matrix elements obey the following symmetry relations 



M-s-,!> = (-1) 



(4.25) 



M 



o 



(a) 



ii 




FIG. 1. Scattering of a plane elastic wave by a plane of spheres: (a): 
incident from the right. 



the wave is incident from the left; (b): the wave is 



V. SCATTERING BY A SLAB 



In what follows we need to evaluate the scattering properties of a slab which by definition consists of a number of 
layers (planes of spheres). For this purpose it is convenient to express the plane waves on the left of a given plane of 
spheres with respect to an origin, A;, on the left of the plane at — d/ from its center, and the plane waves on the right 
of this plane with respect to an origin, A r , on the right of the plane at d r from its center, i.e. a wave on the left of 
the plane will be written as ^2 gi u s si exp [iK g „ • (r — A;)] and a wave on the right of the plane will be written as 

Sgi M gi ex P l^gv ' ( r — ^- r )] ^he relationships between the amplitudes of the incident and of the reflected and 
transmitted waves, when these are expressed with respect to the above origins, follow directly from the corresponding 
equations of Section IV. Accordingly, the amplitudes of these waves are related through the Q-matrix elements given 
below 



Qgijgv = M sks'V exp 
3gi;gV = M+r sH , exp 



gi;g'i> 



exp 



Q 



gY;gV — ^fgijg'i' eX P 



i(K+ •d r + K+ v ,-d i ) 
i(K+ • d r - K gV ■ d r ) 
-i(K~ • d ; - K+„, • d,) 
-i(K" • d, + K gV , • d r ) 



(5.1) 



whose physical meaning is made obvious by their one-to-one correspondence with the M-matrix elements of Section IV 
From this point on, we shall write the above matrices in compact form as: Q 1 , Q 11 , Q m and Q IV . 
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Q'(N) 

— ► 



Q'(N+1) 
► 



Q'(N,N+1) 



Q In (N) 



N))| 



Q ln (N+l) 



^1 



Q In (N,N+l) 



I(q 



(N) 



I 

ICQ"(n+1) 



'i 



|fQ"(N,N+l) 



Q ,V (N) 



Q IV (N+1) 



Q ,V (N,N+1) 



FIG. 2. The Q matrices for two successive layers are obtained from those of the individual layers. 

We obtain the transmission and reflection matrices for a pair of two successive layers, N and N + 1, to be denoted 
by Q(N, N + 1), by combining the matrices Q(iV) and Q(N + 1) of the two layers, as shown schematically in Fig. ||. 
One can easily show that 

Q l (N, N + 1) = Q\N + 1) [I - Q u {N)Q m (N + 1)] ~ l Q\N) 
Q n (N, N+l) = Q U {N + 1) + Q\N + l)Q n (A0 

x [I - Q m {N + l)Q n (A0] ~* Q IV (iV + 1) 
Q ni {N, N + 1) = Q m (A0 + Q IV (7V)Q m (7V + 1) 

x [I - Q n (A0Q m (jV + 1)] _1 Q l (N) 
Q W (N, N + 1) = Q IV (A0 [I - Q U1 (N + l)Q n (iV)] ^ Q IV (N + 1) . (5.2) 
For example, knowing that 

[I - Q n (iV)Q m (iV + 1)] _1 = I + Q n (A0Q m (7V + 1) 

+Q n (iV)Q m (7V + l)Q n (iV)Q m (7V +!) + •••, (5.3) 



we can write the first of Eqs. ( |5.2| ) as follows 

Q r (iV, N+1) = Q I (N+ 1)Q\N) + Q\N + l)Q n (7V)Q m (^ + l)Q r (iV) 

+Q I (iV + l)Q n (A0Q m (Ar + l)Q II (7V)Q m (A r + l)Q r (iV) + • • ■ . (5.4) 

The meaning of the terms is obvious. The first term signifies transmission through the A^th layer, followed by 
transmission through the (N + l)th layer. The second term signifies transmission through the Nth layer, followed 
by reflection by the (N + l)th layer, followed by reflection by the ^Vth layer, followed by transmission through the 
(N + l)th layer. The third and higher terms can be interpreted in the same way: a wave incident from the left 
on the pair of layers will be multiply reflected, any number of times, between the layers before e xiti ng the pair by 



transmission through the second layer. In similar fashion one can understand the remaining Eqs. (5.2). All matrices 
refer of course to the same u> and k|| . We remember that the waves on the left (right) of the pair of layers are referred 
to an origin at -di(N) (+d r (N + 1)) from the center of the Nth {{N + l)th) layer. The choice of d t (N) and d r (N) is 
to some degree arbitrary, but it must be such that A r (N) coincides with A;(A^+ 1), in accordance with the definition 
of these quantities (see Fig. ||) . 
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n rf .\, O A ' (N+1) 



* A,(N) = A,(N+1) 




<jx;n+i) 

d,(N) 



FIG. 3. Putting together a pair of planes of spheres. 

It is obvious that by the same process we can obtain the transmission and reflection matrices of three layers, 
by combining those of the pair of layers with those of the third layer; and that we can in similar fashion obtain 
the transmission and reflection matrices for a slab consisting of any finite number of layers. In particular, having 
calculated the Q-matrix elements of a single layer, we can obtain those of a slab of N max = 2 M identical layers 
by a doubling-layer scheme as follows: we calculate the Q-matrix elements of two consecutive layers in the manner 
described above, then using as units the Q-matrix elements of a pair of layers, we obtain those of four consecutive 
layers, and in this way, by doubling the number of layers at each stage we obtain the Q-matrix elements of the slab. 

In summary, for a plane wave of polarization v' ', £\ [win]^ exp iKg V , • (r — A^) ej, 
incident on the slab from the left, we finally obtain a reflected wave 
J2 S i Kf]~jexp [iK~„ • (r - A L )] e, on the left of the slab and a transmitted wave J2 S i Nr]^ exppK^ • (r - A^)]^ 
on the right of the slab, where A^ (Ar) is the appropriate origin on the left (right) of the slab. We have 



M gi = EC'i'Wg'i' (5-6) 

i' 

where the Q-matrix elements are those of the slab. In the present formulation of the problem we assume that the 
host material between the spheres extends to the left and right of the slab to infinity. However, the extension of the 
formalism to deal with different materials on the left and right sides of the slab can be easily effected by treating the 
interfaces as scattering elements described by appropriate Q-matrices, as in the case of photonic crystals.EZI 

A transmitted beam (a plane wave with a real K^ uz component of the corresponding wavevector) carries with it 
an energy flux density which, averaged over a time period T — 2it/u;, gives (a formal proof of this formula can be 
obtained by applying the standard definition of the Poynting vector P for elastic waves, Pi — —UikUk, to a plane 
waveEa) 




We recall that for a longitudinal wave (y = I) i = 1, while for a transverse wave {y = t) i = 2, 3; and * denotes, as 
usual, complex conjugation. We note that the quantity in braces in Eq. ([T7]) gives the square of the amplitude of the 
displacement associated with the given plane wave. The transmitted energy per unit area of the slab per unit time, 
associated with the g, v beam, is given by the magnitude of the z component, |-Pg^ z | , of Pg V ■ A similar formula gives 
the energy flux associated with any of the propagating reflected beams, or with the incident wave. For the reflected 
beams we have 

Pgi = \p"<? v \ E M gi (W] gi ) * I K" (5.8) 
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And for the incident wave (of given k|| + g' and polarized along the i' direction) we obtain 

ip^,{[ Uin ]+ v ([ Uin ]+ v )*}K+ I/ , 



■pin _ 



(5.9) 



The reflected energy per unit area of the slab per unit time associated with the g, v reflected beam is given by the 
magnitude of the z component, \Pg VZ | , of Pj^ and the incident energy per unit area of the slab per unit time is given 
by the magnitude of the z component, |-P~v z |, of PSj/- By definition the reflectance is given by 



ft(w,k|| +g',0 

and the transmittance by 

T(w,k||+gV) 



y 


prf 1 




pin 
S'i'z 





Egi^M* ([Wrf] gi J K+ 



Jgj 



Jg, 



y 


ptr 1 
guz \ 




r>in 
g'i'z 





(5.10) 



Jg'i' 



Jg'j' 



(5.11) 



where we have denoted explicitly the dependence of these coefficients on the incident parameters. As long as there 
are no energy losses in the slab, we have 



T+n = i. 



(5.12) 



VI. THE COMPLEX BAND STRUCTURE 



We view the infinite crystal as a s equence of identical layers parallel to the xy plane, extending over all space (from 
2 — > — oo to z — > +oo). If Eq. (4.1) is the 2D space lattice for the layer, and a3 is a vector which takes us from a 
point in the iVth layer to an equivalent point in the (N + l)th layer, then {a 1; a 2 , a 3 } is a set of primitive vectors for 
the crystal. 

In the region between the Nth. and the (N + l)th layers the wavefield, of given u> and kn , has the form 

U M = E i u U N ) ex P t iK ^ ' ( r - M^))] + exp [iK~ • (r - A r (N))} } e, . (6.1) 

The coefficients u si (N) are related to the u^(N + 1) coefficients through the scattering properties of the Nth layer. 
We have 



g'i' g'i' 

U + (7V + 1) = E<4;gV<vW + E^;gvV,'(^ + 1) (6.2) 

g'i' g'i' 

where Q are the transmission/reflection matrices of the layer. 
A generalized Bloch wave, by definition, has the property 

u s si (N + 1) = exp (ik • a 3 ) u s gi (N) , (6.3) 

with 

k=(k ] |,A!,(w J k||)) (6.4) 

where k z is, for a given k||, a function of w, to be determined. 

We choose the reduced k zone of reciprocal space as follows: (k||, k z ) where k|| = (k x , k y ) extends over the SBZ of 
the given crystallographic plane, and — |b 3 | /2 < k z < |b 3 | /2, where b 3 = 27ra! x a^/ai • (a 2 x a 3 ) = e z 27r/a 3z . The 
periodicity of the frequency band structure parallel to the xy plane follows from Eq. ( |6.l[ ); for replacing k|| by k|| +g in 
this equation renames the coefficients without changing the form of the wavefunction. Also, because the eigenvalues 
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of Eg. ( |6.5D below are of the form exp (ik - as), values of k z differing by an integral multiple of |lc>3 1 correspond to 
the sam e Bloch wave ; which establishes the periodicity of the band structure normal to the xy plane. Substituting 
Eq. ( |6.3D into Eq. ([T^) we obtain, after some algebra (see Appendix |c]) , the following system of equations 

Q 1 ^ ( U+ ( N ) \- n, \( U+ ( N ) \ (a o 

-[Q^-'Q^Q 1 [Q 1 ^ 1 [I - Q m Q n ] ){u-(N + 1)) - exp ^ lk - a3 ^u-(7V + l) ) W 

where are column matrices with elements: u^ ±1 , u^ l2 , u^ l3 , u^ 2l , u^ 22 , Mg a3 , Ug s i> • • •• In practice we keep 
<?max g vectors (those of the smallest magnitude) in which case are column matrices with 3g max elements. The 
enumeration of the g vectors implied in the above definition of is of course the same with the one implied in 
relation to the Q matrices, each of which has 3(7 max x 3 g max elements; I is the 3<? max x 3g max unit matrix. For given 
kn and u>, we obtain 6<? max eigenvalues of k z from the eigenvalues of the 6<? max x 6g max matrix on the left-hand side 



of Eq.(6.5). The eigenvalues k z (cj;\au), looked upon as functions of real u, define, for each ky, 6g n 
complex k z space. Taken together they constitute the complex band structure of the infinite crystal associated with 
the given crystallographic plane. A line of given k|| may be real (in the sense that k z is real) over certain frequency 
regions, and be complex (in the sense that k z is complex) for oj outside these regions. It turns out that for g iven 
k|| and uj, out of the 6g max eigenvalues of fc z (w;kii) none or, at best, a few are real; the eigensolutions of Eq.([T^) 
corresponding to them, represent propagating modes of the elastic field in the given infinite crystal. The remaining 
eigenvalues of fc z (w;k||) are complex and the corresponding eigensolutions represent evanescent waves. These have an 
amplitude which increases exponentially in the positive or negative z direction and, unlike the propagating waves, do 
not exist as physical entities in the infinite crystal. However, they are an essential part of the physical solutions of 
the elastic field in a semi-infinite crystal (extending from z — > — oo to z = 0) or in a slab of finite thickness. A region 
of frequency where propagating waves do not exist for given k|| constitutes a frequency gap of the elastic field for the 
given k|| . If over a frequency region no propagating wave exists whatever the value of k|| , then this region constitutes 
an (absolute) frequency gap. 

Finally, it is worth noting that, when there is a plan e of mirror symmetry associated with the surface under 
consideration, the eigensolutions (Bloch waves) of Eq.(3.5) appear in pairs: fc z (aj;k||) and — fc z (w;k||). 



VII. AN EXAMPLE 



We demonstrate the applicability of our method by applying it to a specific example, which has also been considered 
by Sprik and Wegdam:E3 a system of silica spheres of radius S = 0.25 pm centered on the sites of an fee lattice with 
a lattice constant of 1 pm; the host material being ice. The relevant parameters are, for silica: p — 2200 Kgm~ 3 , 
q = 5970 ms" 1 , c t = 3760 ms" 1 , and for ice (at -16 °C): p = 940 Kgm" 3 , q = 3830 rns -1 , c t = 1840 ms" 1 . We 
view the crystal as a succcession of planes of spheres parallel to the (001) direction of the fee lattice. Fig. |] shows the 
frequency band structure normal to the (001) plane (k|| = 0) and the corresponding transmission spectrum for both 
longitudinal and transverse waves incident normally on a slab of the above crystal consisting of 16 layers. _ . 

To begin with, we compare the band structure normal to the (001) plane with the results of Sprik and Wegdam £3 
obtained by the plane-wave method (using 343 plane waves) with an accuracy, as stated by the above authors, of a 
few percent. Our results obtained with an angular momentum cutoff irmm = 4 and 13 g vectors are converged within 
an accuracy of 10 -3 , and they agree with those of Sprik and Wegdam E_2l within the stated accuracy of their results. 
Furthermore, the evaluation of the transmission coefficient, easily obtainable by our method but not possible by the 
plane-wave method, confirms the validity of the band-structure calculation. The transmittance curves are shown in 
Fig. U by the shaded curve for the transverse waves and by the solid line for the longitudinal waves. The oscillations 
in the transmittance curves are due to multiple reflections at the edges of the slab (Fabry-Perot like oscillations). As 
expected, for frequencies within a gap the corresponding transmission coefficient vanishes. 

The eigenmodes of the phononic crystal are strictly speaking always hybrid, having a longitudinal and a transverse 
component. However along symmetry directions the following situation may arise. We consider for simplicity the 
normal modes corresponding to k|| = (the direction normal to the surface). In this case the component of the field 
associated with the g = beam is either longitudinal or transverse; the field associated with the g ^ components 
need not be of the same type. However, only the g = component couples to the external field (incident, reflected 
and transmitted waves), if (k|| + g) > for g ^ 0. Therefore an incident longitudinal or transverse wave will excite 
a mode (or modes) in the interior of the crystal with a g = component of the same type. As long as the amplitude 
of the g = component of these modes is much greater than those of the g ^ components, which is the case in the 
example we have considered, the transmitted and reflected waves will be of the same type as the incident wave, but 
this need not be the case in general. Our results shown in Fig. are termed longitudinal or transverse in the above 
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sense. For waves incident at an angle on the surface of the slab (k|| ^ 0) the above distinction between longitudinal 
and transverse waves no longer applies (the g = component of the elastic field inside the crystal is a hybrid one) 
and therefore an incident wave of a specific type (longitudinal or transverse) will give rise to reflected and transmitted 
waves of a mixed type. 




-l o no 0.5 l 

kal2n Transmittance 
FIG. 4. The phononic band structure at the center of the SBZ of a (001) surface of an fee crystal of silica spheres in ice 
(a); and the corresponding transmittance curve of a slab of 16 layers parallel to the same surface (b). The lattice constant is 
1 fira and the radius of the spheres is 0.25 (im. In (a) the black lines represent longitudinal modes (in the sense defined in the 
text), the grey lines transverse modes, and the dotted lines are deaf bands. Correspondingly in (b) the solid line shows the 
transmittance for longitudinal incident elastic waves; and the shaded curve that for transverse incident waves. 




r x M r 

FIG. 5. Projection of the frequency band structure on the SBZ of the (001) surface of the fee phononic crystal described in 
the caption of Fig. ^. The shaded areas show the frequency gaps in the considered frequency region. The inset shows the SBZ 
of the (001) surface. 

In Fig. H we show the projection of the frequency band structure on the SBZ of the (001) plane along its symmetry 
lines. This is obtained, for a given ky , as follows: the regions of lu for which there are no propagating states in 
the infinite crystal (the corresponding values of all k z (w,k||) are complex), are shown shaded, against the white 
areas which correspond to regions over which propagating states do exist (for a given u there is at least one solution 
corresponding to k z (£j,k||) real). We note the existence of a narrow absolute gap, denoted by Awg, extending from 
2.82 GHz to 2.89 GHz. An absolute gap at approximately the same frequency and of approximately the same width 
was found by Sprik and Wegdam.113 

A considerable number of bands of longitudinal and transverse waves exist above the absolute gap (see Fig. ||). 
Below this gap we have two bands of transverse waves, which are doubly degenerate, extending from GHz to 1.77 GHz 
and from 2.01 GHz to 2.82 GHz, with a gap in between. On the other hand, a non-degenerate band of longitudinal 
waves extends from GHz to 2.48 GHz. In addition to these bands, we find a non-degenerate band, extending from 
2.40 GHz to 2.82 GHz, for which the g = component of the corresponding eigenmodes of the elastic field vanishes. 
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Because the g = beam is the only one which matches (couples with) a propagating wave outside the crystal, an 
internal mode with a vanishing g — component is not excited by the incident wave. Therefore, if this were the only 
band over the stated frequency region, the wave would be totally reflected. However, in our example transmission 
through the slab in the frequency range of this deaf band occurs, because other bands with non-vanishing g = 
components exist in the same frequency region. We note that analogous deaf bands are known to exist in photonic 
crystals .Ej 

The long wavelength limit (k z — > 0) is represented by the linear segments of the dispersion curves, the slopes of 
which determine the propagation velocities of longitudinal and transverse waves (cj = 3893 ms" 1 , Ct = 2033 ms -1 ) in 
a corresponding effective medium. 



VIII. CONCLUSION 



We have shown that for a system of non-overlapping elastic spheres arranged periodically in a host medium of 
different elastic coefficients one can, using the formalism of the present paper, calculate accurately and efficiently the 
phonon spectrum of the infinite crystal and, also, the transmission, reflection, and absorption coefficients of elastic 
waves incident on a slab of the material of finite thickness. 



APPENDIX A: 



The nonzero elements of the T matrix for a solid sphere in a solid host are 



T 

J- 1 



MM 



T, 



NN 



(p s zf I 'px 2 )je(z t )[x t j' e (x t ) -jejxt)] - je(xt)[ztj'e(z t ) - je(z t )] 
ji(x t )[zthj (zt) - hj(zt)] - (p s zf I px 2 )hj (z t )[xtf e (x t ) - je(x t )} 
W e NN 



-5ffi5 n 



e > i, 



-Su>5n 



£,£' > 1, 



NL 



{zt/zi)— ^ o u ,5 mm , , 



T, 



LN 



(zi/zt) 



e> i, £' > o, 

, I > 0, I' > 1 , 



J- 1 



LL 



LL 



-Sei'Smm' , £,£' >0 , 



(Al) 



with 



Sq v , x v — Sq sv and v = l,t. The elements of the 4x4 determinant Dt are given by 



(hi 

d 3 i 
(hi 
d\2 
d 2 2 

dz2 

di2 



z t hf (zt) + hf(zt) 
£(£+l)h+(zt) 

[£(£ + 1) - z 2 t /2 - 1] h+(z t ) - z t h+\z t ) 

£{£+!) [zthj'(z t )-hj(z t ) 

ht(zi) 
zihj (zi) 

zihj (zi) - hf(zi) 

[£{£ + 1) - z 2 t /2] hf(zi) - 2z l hf(z l ) 



(A2) 



du = x t j' t {x t ) +je(x t ) 



d 2 3 

d33 
<^43 

du 
d 2 4 
d 3 4 

C?44 



: { Ps zllpx\) { [£(£ + 1) - x?/2 - 1] j t {x t ) - xtjfat)} 
: {p s z 2 /px 2 t )£(£ + 1) [x t f e {x t ) - j t (xt)] 
-ji(xi) 
■ xij' e (xi) 

: (p s z 2 /pxf) [xif e (xi) - je(xi)} 

■■ (p s z 2 /px 2 ) { [£{£ + 1) - x 2 /2] j t (xi) - 2x a 'M)} , 



where j' e and h\ denote the first derivatives of the spherical Bessel and Hankel functions, respectively. Wf p are 
given by the following determinants 
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< 




di3 


du 




d 2 2 


d 2 3 


d 2 i 




d 3 2 


d33 


d 3 4 




di2 


di3 


da 



NL 



d\ G?12 d\3 du 

d 2 d 2 2 C?23 <^24 

^3 d3 2 ^33 C?34 

d\ di 2 dn3 da 



W[ N = 



du d 1 di3 di4 

dii d 2 d 2 3 d2A 

d 3 l d 3 C?33 C?34 

c?4i d^ (I43 da 



TT', 



LL 



du d\ d\3 c?i4 

d 2 i ^23 ^24 

^31 dg (i33 C?34 

£^41 (^4 C?43 C?44 



(A3) 



where 



rf i = z tf e ( z t) +ji(zt) 
<ff = £{£+!) j t {z t ) 



-,N 
J 3 
,N 



[£(£ + 1) - - 1] Mz t ) - ztj'^zt) 



(A4) 



d% =£(£+l)[z t j' e (z t )-j e (z t )} , 



and 



di =ji(zi) 

d 2 = z d'e( z i) 

^3 = z d'i( z i) -Je(zi) 

d\ = [£(£ + 1) - zf/2] u{zi) - 2zif i (z l ) 



(A5) 



APPENDIX B: 



A longitudinal (transverse) spherical wave about R„ ^ remains a longitudinal (transverse) wa ve w hen expanded 
about the origin of coordinates (R n = 0), and therefore the matrix elements of l~i defined by Eq. ( |4.9| ) are obtained 
independently for longitudinal and transverse waves. 

For the transverse waves the evaluation of these elements proceeds as in the case of the electromagnetic (EM) field 
described in Ref. 37. We note that the M transverse elastic wave corresponds to the H component of the electric 



field of the EM wave and the N transverse elastic wave corre spon ds to the E component of the electric field of the 
EM wave (compare Eqs. (15) and (16) of Ref. ^ with Eqs. (4.7) and (4.8) of the present article). Therefore, the 
fl pp ' -matrix elements with P = M,N and P' = M,N can be taken directly from Re f. |37| . Taking into account the 
fact that the expansion coefficients in Eq. (ETjj) above and those in Eq. (15) of Ref. |37| are multiplied by different 
constants, one readily obtains 



MM 
£m:£'m' 



{£{£+!)£' {£' + !)} 



1/2 



r^NN 

iL £m;£'m' 



QMM p pi > 1 



£J' > 1, 



(Bl) 



where 



n MN 



-n 



NM 



-1/2 



= {2£+\)[£{£+\)£'{£ l + 
x|(87r/3)V 2 (-l)™ a y , 4-^+ 1 (ft)^-i,r n +i(l > -l^m) 

+m'Z e e ^r(q t )[(£ + m)(£~m)/(2£-l)(2£+l)] 1 / 2 }, £,£' > 1, 



?im'(<lt)= exp(ik|| ■R n )Gt m .if m >{-R n ;qt). 

R„^0 



(B2) 



(B3) 
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G Wm ^-R ri ; gt ) = 47r^(-l)^'- r )/ 2 (-l) m '+ m ''^ m (/m';rm' / ) 



xh+(q t R n )Y e 7 m (-R„) : 



(B4) 



B <ra (ffn';fm")E / df Y t m (r)Y e T (r) Y^ m (f ) . 



(B5) 



The expression for £l MN can be simplified further by the evaluation of the Bg m coefficients defined by Eq. (B5). Using 
standard formulas (see, e.g, Ref. |38|) one finally obtains 



o — m' rn ryt— lm- 1 / \ , _/Amv^— lm / \ , o rn —mryl—lm+lr \ 

-2a e , Tf^'m'-i ( gt ) + mC £ m ^£w (<7t) + 2a£? 7^ £</ m / + i (?t) 



^m-f'm' — (2^+1) 



1/2 



atm _ _ q mn p pi > 1 



(B6) 



where 



Tf = ^ + m)(i + m - l)] 1/2 /[(2^ - 1)(2€ + 1)] 1/2 and 



C - [(i + m)(£ - m)] 1 / 2 /[(2^ - 1)(2* + 1)] 1/2 . 
The derivation of the above formulae is based on the following relationEll 

h+(qr n )Y e m (r n ) = ^ Gv B '(-R»;#(jr)17'(f) ; 



(B7) 



(B8) 



which expresses a scalar spherical wave about R n 7^ 0, as a sum of spherical waves about the origin (R„ = 0). The 
longitudinal wave, described by the third term of Eq. (4.7), is obtained by multiplying Eq. (B8) with exp(ik|| • R„) 
and summing over all R„ ^ 0, which immediately tells us that 



(B9) 



The evaluation of the matrices 17 involves the evaluation of the matrix Z which is a well known quantity in the 
theory of low-energy electron diffraction (LEED) and a computer program for its evaluation is already available in 
the literature!^ Further calculation is made simpler by the following property of Zf™ (<?„)Ej 



Z(™'(q u ) = 0, unless £ + m + £' + m! : even. 



It follows from Eq. ( |BTo| ) that 



Sl%$, m , = SlZZtw = Vim-A'm' = 0, unless i + m + e + m': even, 



AIN 



n 



NM 



0, unless £ + rn 



m : odd. 



(BIO) 



(Bll) 



APPENDIX C: 



Eq. (3.5), initially derived by McRae for th e ca se of electron scatteri ng b y atomic layers can be proven as follows. 
We replace the quantities on the left of Eqs. (6.2) with the aid of Eq. (6.3) to obtain 



explik ■ a.ji ( *„ %• 



u+(N) 
u-(N + l) 



Q 1 Q n 
I 



u+(N) 
u~(JV + l) 



The inverse of the matrix on the left of Eq. (CI) is given by 



I 



Multiplying both sides of Eq. 



-[Q IV PQ m [Q IV ]" 
with the above matrix we obtain Eq. (|6.5|). 



(CI) 



(C2) 
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